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Abstract 



We study the tautomeric transitions in base pairs of DNA considering clastic 
properties of DNA as classical and tunneling of protons as quantum, and show 
that the dynamics of the transitions admits of soliton like solutions whose shape 
and size strongly depend on the structure of the double helix. In particular, we 
have found that the set of discrete breathers can be drastically modified by 
the interplay of the torsional and clastic constants. Our results may have a 
bearing upon substitution mutagenesis within the framework of Watson-Crick's 
approach, and in this respect the breather soliton could describe conformations 
corresponding to point mutations. The numerical simulation of soliton dynam- 
ics suggests that an initial distribution of base pairs with low probability of 
mutation per pair but of a sufficiently large number of base pairs involved, 
could move and gather around a site so as to form a set of base pairs with high 
probability of mutation, for a period of time approximately 1 fxsec. We suggest 
that the irradiation of DNA at frequencies of the proton tunneling, that is in 
infra-red region, could cause mutations. 
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1. Introduction 



According to the Watson-Crick hypothesis, (!]], the life is built around a 
symmetrical figure, the double helix of the DNA molecule (see Fig.l), which 
comprises the two strands linked together by purine-pyrimidine base-pairs of 
adenine-thymine (AT) and guanine-cytosine (GC), the four chemicals A,T,G,C 
existing in various isomeric forms, or tautomers, that may change into one an- 
other (see Fig. 2). Under ordinary conditions the equilibrium shifts towards the 
amino-form for adenine and guanine, and the keto-form for thymine and cyto- 
sine. But the imino-form for the adenine and cytosine, and the enol-form for 
guanine and thymine are also possible, even though rare; in fact, they corre- 
spond to concentrations of to 10~ 5 moles/liter, Q. To see the implications 
wrought by the tautomeric transitions let us recall that the sequence of base- 
pairs constitutes the genetic information of cell. It should be noted that adenine 
will pair only with thymine by two hydrogen bonds, and guanine with cytosine 
with three hydrogen bonds, so that exact copies of the DNA are produced dur- 
ing the replication (see Fig. 3). But the complimentarity between the bases is 
completely changed if a tautomeric transition takes place; in fact, owing to the 
structure of hydrogen bonds other combinations become possible (see Fig. 4) , 

Aimino * * G , A < ► Cimmo (1) 

Genol * * T : G < ► T eno i 

in contrast to the usual and stable ones 

A< >T , G< >C 

Another opportunity for generating " unnatural" pairs arises from the tunneling 
of protons in hydrogen bonds (see Fig. 5), which results in the formation of the 
pairs 

{A < > T) ^ (^-immo < * T^enol) (2) 

(G < ► C) )' {Genol * * Ciraino) 

During the replication, tautomeric transition driven by the proton tunneling in 
conjunction with the complimentarity according to (j2j) may lead to the change 
of base-pairs 

(A<— T) =► (G<— >C) (3) 
(G < ► C) =► (A < ► T) 

and result in loss, or corruption, of genetic information, i.e. mutations, ||, Q. 
The specific case given by the diagram (|||) is called transition mutations; it has 
the property of being reversible, i.e. able to go back to the wildlife type. 

The arguments given above constitute the main points of the theory of spon- 
taneous mutations suggested by Crick and Watson, ||, ||, Q. It is based on 
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the assumption that the transitory tautomeric shifts of base-pairs may occur 
during the replication, i.e. when two molecules of DNA are formed from a 
paired molecule, so that the double-stranded molecule is split into two single 
strands, each of which controls the synthesis of a new strand complimentary 
to itself with the help of the special enzyme called DNA polymerase. It has 
been realized that the latter plays an active role in the selection of bases at 
replication, 0, so that it may affect the mutation rates. Thus, tautomeric tran- 
sitions are not a unique cause of mutation; the situation is more subtle, and 
many questions, of quite a classical nature, wait their solutions. Nonetheless, 
the original idea of Watson and Crick still conserves its appeal, and even more 
so as its new links with other phenomena related to the mutagenesis are brought 
to light (see ||). So, Robinson et al, || report that the enol tautomer of iG, 
that is 2'— deoxyisogine, may form at physiological temperature (37°) and pair 
with thymine in a Watson-Crick geometry (see Fig. 14); thus, iG present as the 
nucleoside, results in the formation of incorrect base-pairs during in vitro repli- 
cation, Robinson et al, §, suggests that iG ■ T pairing 
may have a bearing on mutagenesis in vivo involving tautomers of the common 
nucleobases. On the other hand, Fresco et al, pjfl , have found that the imino 
tautomer H0 5 dCyt may serve as an example of an unfavored base tautomer 
making for substituting mutagenesis. 

Mutations within the framework of the Crick- Watson model of DNA and 
in conjunction with the concept of tautomeric transition, have been drawing 
attention, beginning from the early fifties, Q, ||, @, to the present time, and 
involved the use of condensed matter theory. So, one of the first papers in 
this direction was published by Geracitano and Persico, who suggested 

that there should be expected a collective behavior of codons, resembling that 
taking place in hydrogen-bonded ferroelectric crystals. 

In this paper we intend to look after the interplay between tautomeric 
transitions in base-pairs and elastic properties of the double helix. Since the 
7T— electrons of the tautomeric rings of the nucleotides have direct bearing on 
the interaction of the plates of adjacent base-pairs, jl7), E3, we suggest that 
the tautomeric transition of base-pairs should substantially influence the dis- 
tribution of delocalized electrons of the nucleotides, i.e. the 7r— electrons, and 
result in deformation of the elastic system of DNA. It is worth noting that tau- 
tomeric transitions may occur in several base pairs, not necessary adjacent, at 
a time, and their dynamics is determined by the proton tunneling. For one 
thing the latter is determined by the electrostatic interaction, i.e. the dipolc 
forces, between the protons belonging to adjacent base-pairs, and for another 
by the elastic system of the DNA molecule, which should play a role like that 
of the crystalline lattice of the polaron theory. The situation is similar to that 
considered in the Davydov theory for the a— helix of proteins, To put 

these arguments in a more quantitative form we begin by recalling certain facts 
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concerning the elastic properties of DNA. 



2. The torsional modes of DNA 

From the mechanical point of view, the DNA molecule is a very unusual 
object. For one thing one can visualize it as an elastic rod and consider it 
within the framework of the mechanics of continuous media, thus obtaining 
reasonable agreement with experiment, ppj , for another it is often necessary to 
study DNA as a discrete system similar to crystalline lattice. In the present 
paper we have to deal with an intermediate situation for which a cautious use of 
the methods of continuous media is justifiable. In fact, as far as the transport 
of torsional stress (torque) along DNA is concerned, its estimates obtained by 
various means diverge widely. The numerical values derived with the help of the 
theory of continuous media, ^2|, are of the order r oc 10 -17 dyne ■ cm, [^0|, pl[ , 
whereas there is the experimental evidence, p3| , that it can attain the value of 
r oc 10~ 13 dyne ■ cm. Philip Nelson, po| i suggested that these deviations could 
be due to small bends in the helix backbone, so that one may assume 

r oc 1CT 17 ~ 1CT 13 dyne ■ cm 

As was mentioned above, the interplay between the torsional stress due to the 
relative motion of the base-pairs and the proton tunneling is very important. 
We shall use the approach worked out in J|4|] and [^(| to describe the elas- 
tic properties of the double helix. Thus, the double helix is considered as a 
one-dimensional lattice of vectors y n describing the mutual position of the two 
strands at sites corresponding to the base-pair of index n. It is important that 
the system has a twisted ground state characterized by the twist vector fi, so 
that the elastic energy of the molecule can be cast in the form 

N N N 

H tor = Y,-M{dmf + Y,\ Ki yy*f + E^' 2 ( 4 ) 

i—1 i—1 i—1 

where the first term is the kinetic energy, the second one the elastic torsional 
energy and the last one corresponds to the separation of the two strands. The 
covariant derivative that accommodates the torsion of the molecule, reads 

W, = - (tfi+i - Vi + tixyij 

Here a is the spacing between the adjacent nucleotides, M is the mass of base- 
pair. It should be noted that we are considering a very simplified model and 
assume that all sites, corresponding to base-pairs are identical. The subtle 
question is the value of the elastic constant K; obviously enough it has a direct 
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bearing on the torque r mentioned above, and therefore its estimate may read 

K oc 1(T 13 10~ 17 erg 

It should be noted that the calculations within the framework of molecular 
dynamics, (see paper [^7) and references therein), give the upper value for K, 
i.e. close to 10 - 12 -j- 10~ 13 erg. For the sake of simplicity, in this paper we 
shall assume that the torsion vector f2 is always parallel to the axis Oz, that is 

Q = (0,0,0) 

and the vectors y n describe only transversal motions, that is y 3 = 0. 

As was mentioned above the tautomeric transitions are driven by the proton 
tunneling, and therefore we shall describe them quantum mechanically, that is 
the stable amino/keto form corresponding to the ground state of proton, and the 
unstable imino/enol one to the excited state, [ p5j . In accord with the qualitative 
character of our approach we neglect the fact that the tautomeric transitions 
in question involve the tunneling of more than one proton, and assign only one 
proton to each site of the lattice. 

It is important that there are few hydrogen bonds in which the protons are 
transferred towards the imino/keto groups, or if one uses the concept of the 
two-level system, the excited states. Therefore, one can consider the system 
as being close to equilibrium, or only weakly excited. This suggestion is very 
important for what follows. 

We shall describe the states of a base-pair at site n with the Bose operators 
b£, b n that verify the usual conditions 

[ b n, &m] = b n b m ~ b m b n = S nm , [b„, b m ] = [&+,&+] = 

so that the energy of the protons, neglecting the interaction with the elastic 
degrees of freedom, reads, (25) 

Hp = J2 E b+b n + nJ2Kb n+ i + b+ +1 b n ) 

n n 

Here E a is the energy of the tautomeric shift; its estimates depend on the choice 
of nucleotide and according to quantum chemistry calculations vary within the 
range of 2 -=- 10 Kcal, (see g] and references therein). The constant k could be 
ascribed to dipole interactions between adjacent sites, similarly to Davydov's 
theory, |19J . Presently, there are no reliable estimates of its value (see below); 
by analogy with the Davydov theory one may assume that it should correspond 
to the characteristic frequency of proton , or tautomeric, excitation of the order 
10 , or less. This figure is generally accepted (see below). 
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The central point of the model introduced in p5[ is the interaction between 
the elastic degrees of freedom of DNA and the tautomeric transitions, or the 
proton tunneling in nucleotides; it reads 



An argument in favor of this choice is that it takes into account the deformation 
of positions of adjacent base-pairs and thus its influence on the tt— electrons of 
the bases, and therefore, the tautomeric transitions, or the related excitations 
of protons. According to the theory of [Q, the interaction could be apprecia- 
ble. Thus, one may suggest that the interaction term could be larger than the 
tunneling term in the equation for Hp given above. 

Concluding we may state that the total energy within the framework of the 
model introduced in [E5| has the form 



Htotal = H tOT + Hp + Hj 



The conditions discussed above are exactly those used for Davydov's theory, 
[ [l9| . Let us recall its main points. It is assumed that the state of the system 
can be described by a trial function that has the form 



\V>=^A n (t)-b+\0> (5) 

n 

where > is the vector designating the ground state of the system, that is all 
the base-pairs, or the protons in the hydrogen bonds, being in the ground state. 
The amplitudes A n {t) are subject to the constraint 

Ei a «wi 2=1 ( 6 ) 



At this point it should be noted the vectors y n describe the dynamics of base- 
pairs, that is relatively massive objects, and therefore one may consider them 
as classical fields, p5| ], [ p6[ . We can derive the size of characteristic frequencies 
for y n from expression (mj of the elastic energy. In fact, the mass M is that of 
the base-pair, that is of the order 500 Dalton, and K is of the same order of 
magnitude as r. Hence, we get the characteristic velocity v for the y modes 



v oc 



G 



Interesting numerical values for the velocity v follow from the equation indicated 
above and the rough estimates for r or K we have mentioned. Indeed, for 
K oc W- 17 dyn ■ cm or less we obtain 

v oc 10 2 cm/sec 

For wavelengths of a few tens of A it gives the characteristic torsion or phonon 
frequencies of the order 

Lo y oc 10 8 ~ 10 9 Hz 

On the other hand, if we use the values for K provided by the molecular dynam- 
ics simulations, |27|| , we get the velocity of excitations of the order lOOOm/sec, 
and u> y oc 10 11 ~ 10 12 Hz, as for ordinary condensed media. It is instructive to 
compare the values of u> y with the transition frequencies for tautomeric reactions 
inside the nucleotides, 

K 

UJp= 2^h 

The estimates for the latter differ considerably, (32|, |33| 

uj p oc 10 6 + 10 11 Hz 

The lowest estimate, 10 6 Hz appears to be not unreasonable (V. Benderskii, and 
J.L.Leroy, personal communications). 

The relative sizes of uip and ui y are important for choosing the right approx- 
imation for the model. In fact, if we are at the lowest end of the spectra cjp, 
then according to the estimate for cu y obtained above the characteristic times 
for the acoustic modes are at least by an order of magnitude smaller than for 
the protons. In this case, we may suggest that the elastic system should fol- 
low the motion of the protons in hydrogen bonds, adjusting itself to it, so that 
the situation is similar to that of the Born-Oppenheimer approximation in the 
atomic theory. 

Thus, we assume, as in paper J25|, that the adiabatic approximation is valid, 
and therefore we may neglect the kinetic energy of the elastic system and take 
into account only its potential energy generated by the field y n . Then we are in 
a position to apply Davydov's method, [^9|, that is to calculate the mean value 

U eff =<V\H tor + H I \V> 

find the minimum, of U e ff with respect to y„, substitute it into the equation 
for the total energy H tota i so as to get the effective Davydov hamiltonian Hd, 
which depends only on the operator variables b n , the classical variables y n 
having disappeared through the minimization. Thus, we obtain an equation 
that has the form of the Schrodinger one 

m^-\V>=H D \V> (7) 
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and in which the wave function \T> > should be of the form (||). The assumption 
that the excited states correspond to the set of two-level systems is accommo- 
dated by the requirement that the operators 6+ are allowed only in the first 
power. It results in a system of equations, called the Davydov equations, for 
the amplitudes A n , which one obtains on equating the coefficients at 6+ on 
both sides of (Q), ( see jl9) for the details ). All the necessary calculation for 
the model of the DNA we employ, had been done in [p5| , on which the present 
paper relies. 

The Davydov hamiltonian for our problem reads 



n n 

n n 

A 2 ea 2 



m. , n 

i2 cn 1 



2k m 2 

ra.n 

and the equation for the amplitudes A 



^cos 1 " 1 "" 1 ^ • cos [(m - n)fo - a)} \A m \ 2 \A n f 



+ E ^ m - n ^ ■ cos [(m - n)(0 - a)} \A n \ 2 b+b. n 



ih— A n = E A n - n{A n+1 + A n -i) 
X 2 A 2 

m 

+ TlJ ( E cos^- m ^-cos[( mi -m 2 )^-a)]\A m A 2 \A m2 \ 2 )A 

mi .7712 

+ y|J(E ^ |m - n| • cos [(m - n)(0 - a)] \A m \ 2 )A n 

rn 

see p5l for the details. 



3. The numerical simulation 

We introduce the reduced variables B n according to the equation 
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A n = e-^ Eot B n {t) 
and cast the equation for A n in the form 



d A 2 
ih—B n = -n{B n+1 + B n - 1 )-—\B n \ 2 B n (8) 



k 

m 

A 2 ea 2 

mi ,m2 



( cos^-^^B^B^Bv 



k kn 2 



Here 



arctan 17 



The non-local terms are a consequence of the structure of double-helix, often 
neglected in considering the dynamics of DNA, pS] . 

Introduce the characteristic frequencies 

" P= 2^' " T= 2^' ^=2^ (9) 
and the dimensionless time 

T = t ■ LJp 

It should be noted that the frequencies u> y and Lo t0 r ar e not identical, io y ^ c<j for . 
Then the Davydov equation takes the form 

,d_ 

-W(J2 \B m \*)B n 



'^rBn = -(B n+1 + B n ^)-W\B n \ 2 B n (10) 



in which 



+WA( E c OS l mi - m2 l0|S mi | 2 |S m2 | 2 )S„ 

mi ,m2 

+^A(^c OS |m ~ n| 0|5 m | 2 )£„ 



W = ^ T (11) 

ea 2 . , 

= Fff (12) 
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Now we aim at making the numerical simulation of equation ( JTo| ) for various 
values of the parameters W, A, looking for solutions of the soliton type. We use 
the term soliton in a sense close to that used by applied scientists, i.e. a solution 
different from zero in a finite region of space, whose size we shall call the size of 
soliton and preserving its shape for very long periods of time. For some values 
of W, A it has the form identical to the usual one, i.e. corresponding to the 
non- linear Schrodinger equation, but generally our solitons are different. The 
standard definition suggests that it be of the form 

Y(x, t) = e l< - qx -^ ip(x - vt) (13) 

in which i\) is a real function. It is by no means clear that solutions that we 
suppose to be solitons, always have the form given by equation ([l3"|). 

The parameter A is a quantitative characteristic that enables us to take 
into account the structure of the double helix, and also the relative size of the 
torsional and deformation energies. In fact, A determines the magnitude of the 
nonlocal terms in equation (|Io|), and in this respect it is worthwhile to note that 
for certain values of A and W we have not been able to find soliton solutions, 
e.g. A = 0.2 and W — 2, at least for physically reasonable sizes of solitons, 
i.e. less than 100 base pairs. But it is important that generally the condition 
A ^ does not forbid the existence of solitons, and its influence only results 
in the size of soliton becoming larger, which is quite natural, for A represents 
non-local terms in equation ([n]). The general case of soliton with A not equal 
to zero, even though small, is illustrated in Fig. 8. 

To illustrate the general situation let us consider the two special cases (for 
the details of calculation see Appendix). 

1. Stationary solutions in the sense that the absolute value, \B n {t)\ does not 



depend on time. For the usual solitons given by (10) this requirement means 
that the velocity v = 0. The typical case is illustrated in Fig. 9, for W — 10 and 
A = 0.5. The half-width of soliton is equal one spacing between base-pairs, that 
is the solution is extremely narrow, and according to our main hypothesis it must 
correspond to the tautomeric transition of a base-pair. The very interesting case 
is illustrated in Fig. 10, W — 5 and A = 0.5. There is a central peak of half-width 
1.5 • a which stands still, and two symmetrical wave packets, moving in opposite 
outward directions. The distance traveled by these wave packets during 0.018 
msec is equal to 33 base pairs. 



2. The usual solitons given by (10). The half-width of these solitons may be 
several tens of base-pair spacings, and thus they could correspond to tautomeric 
transitions taking place in adjacent base-pairs. The typical cases are illustrated 
in Fig. 11, 12. The distance traveled by soliton in Fig. 11 during 0.53 msec is 
equal to 707 base pairs, and the distance traveled by soliton in Fig. 12 during 
0.577 msec is equal to 491 base pairs. It is interesting to note that these solitons 
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move, even though slowly. Hence, one might suggest the picture of tautomeric 
transitions moving along the DNA-molecule. 

Both types of solutions indicated above are stable with respect to perturba- 
tion of W and A. 

Perhaps, the most characteristic feature of discrete non- linear Schrodinger 
equation is solutions that periodically oscillate in time and decay exponentially 
in space, or breathers, |3l| . From a purely qualitative point of view the existence 
of breathers can be inferred from a truncated version of equation (|Io| ) . Let us 
neglect all the terms on its RHS except the first two, that is consider 

ih-^T = -(B n+1 + S n _ x ) - W\B n \ 2 B n 
at 

and look for B n such that 

Dri ^ llji 

a n being real. Next, cast the equation for a n in the form 

-(a„+i - 2a„ + a„_i) - [Wa\ + (2 - e)]a„ = 

Suppose that the soliton we are looking for is large enough so that we may 
change the expression a„+i — 2a n + a„_i for the second derivative. Thus we 
obtain the equation 

a" + [Wa 3 + (2 - e)]a = 
or the conservation law for one dimensional motion with the effective potential 

The soliton solution exists for e > 2, and its size tends to infinity as e — > 2. On 
the other hand for large W we may expect thin solitons. 

The key point is that the nonlocal terms generated by the double helix bring 
serious modifications to the picture given above. 
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We may infer from the examples given above that the dimensionless con- 
stants W and A play a crucial role in determining the form of solitons for equa- 
tion (|Toj) . The general situation to the effect is illustrated in Figs. 16 (a)-(b), in 
which the horizontal axis corresponds to values of is, that is the soliton frequency 
measured in units of uip. It should be noted that v is well defined for solitons 
of the form given by equation ([l3|); in contrast, there is a fine structure in the 
frequency spectrum of breathers, [|l], so that v turns out to be only a rough 
characteristic. In fact, in Figs. 16 (a)-(b) the values of v are determined to within 
about one hundredth of top. As is shown in Figs. 16 (a)-(c) for A = 0, 0.1, 0.15, 
respectfully, the set of W and is, for which there are solitons or breathers, con- 
sists of a line corresponding to breathers, and a region, or domain, for moving 
solitons. The line serves also as a right border for the region of moving solitons. 
The lower border of the soliton region is not strictly defined owing to the fact 
that there are solitons for values of W and is lower than the borders but of sizes 
greater than 100 base pairs, that is outside the physical context of our problem. 
The upper left part of the border is determined by solitons turning out to be 
unstable for values of W and is beyond the boundary. We see that the soliton 
region is decreasing as A grows, and for A = 0.2, Fig. 16 (d), there are only 
breathers, at least under the constraint of their size being less than 100 base 



pairs. It is worth noting that the equation (10) derived in |25j is valid only for 
small A. 

We would like to draw attention to a class of solutions that are not solitons, 
but nonetheless may have a bearing upon the dynamics of tautomeric transi- 
tions. A solution of the type is illustrated in fig. 13. It is characterized by an 
initial set of amplitudes B n (t) which is a broad distribution of the size of 80 
base-pair spacings; after the period of time 0.017 msec, it focuses itself on a 
narrow peak of half-width of one spacing. The peak exists for the brief period 
of time 0.002 msec, and next breaks down into a broad distribution again, i.e. 
a kind of partial self focusing is taking place. Thus, there may exist low prob- 
ability tautomeric transitions distributed over wide areas of the molecule, and 
which may collapse into a small region of the molecule, and stay there for a 
period of time, brief but perhaps sufficient to cause mutation. 

Finally, we wish to tell that our simulations used the standard numerical 
methods, i.e. the trapezoid, the Adams-Boshoft and the Adams-Moulton of the 
fourth order, the algorithm for stiff systems. An important test has been the 
conservation of the normalization condition 

Ei b ««i 2 = 1 ( 14 ) 

n 

For testing the precision of our algorithms we have also used calculations back- 
wards in time. 
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4. Conclusions 

As was shown above the dynamics of tautomeric transitions in DNA depend 
on elastic properties of the latter and proton tunneling in base pairs; W and 
A serving as indicators for possible regimes. Our numerical simulation suggests 
that the interesting tautomeric dynamics may happen for W > 1. This allows 
for sufficiently wide range of material constants of DNA so as to hope the 
phenomenon's taking place. The second constant, A, provides a quantitative 
characteristic for the part played by the double helix; it can totally modify the 
structure of solitons corresponding to tautomeric transitions. 

Depending on the value of W one may expect the existence of two quite 
different kinds of soliton dynamics. The first one belongs to solitons that move at 
velocity of several 0.01 cm/ sec, have a size of several tens of base-pairs, and the 
second of stationary solutions, or breathers, that have a form of peaks over a few 
base-pairs. We may suggest that the second type of solutions correspond to point 
mutations, whereas the first one may describe tautomeric transition moving 
along the chain of double helix. According to the Crick- Watson approach, there 
may happen mutations related to the transition. Thus, one may suggest that an 
action imposed on a set of nucleotide in a region of the molecule might generate 
mutations in a different region owing to the motion of excitations corresponding 
to proton tunneling. 

It is alleged to be known that by substituting the "artificial" nucleotides 
instead of the natural ones, e.g. brom-uracil for thymine (see Fig. 6, 7), one can 
increase dramatically the rate of mutations; this could be due to the increase of 
tautomeric transitions inside base-pairs. At any rate, it is worthwhile to study 
the interplay between the rate of such transitions and mutations. Within the 
context of the present paper, artificial DNA of this kind could ease the stringent 
constraints imposed on W, as was indicated above. 

It is worth noting that the "focusing" of solutions (see Fig. 13) may have a 
very important bearing on mutations. In fact, it amounts to the possibility of 
a weak external influence generating a low amplitude distribution of mutation 
sites that would focus itself later on a high amplitude distribution concentrated 
in a different region of the molecule. Thus, one may expect generating mutations 
by low intensity agents distributed in a region of the molecule, or to put it the 
other way round, acting on a set of codons different from those that suffer the 
actual mutation. 

Finally we would like to point out that the irradiation with electro-magnetic 
waves at frequencies of the proton tunneling in hydrogen bonds may result in 
tautomeric transitions of nucleotides, and, according to the arguments given 
in this paper, mutations. This circumstance could be used for experimental 
verification of our model. 
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Appendix - Radiation Filtering (RF) Algorithm 

The method used in this paper for constructing soliton solutions is based on 
the following heuristic argument. Suppose we have a solution that is close to a 
soliton, then we may visualize it as a central peak that radiates waves of small 
amplitude (see Fig. 10), during the evolution in time. In case of soliton proper, 
the radiation is absent, so that the radiation is a specific feature that gives a 
measure of not being a soliton, so that the difference between the central peak 
and the soliton is carried away by the radiation. Hence, the soliton solution may 
be obtained by annihilating the radiation through the use of an appropriate filter 
(see Fig. 14). 

The actual algorithm runs as follows: Take the parameter of cutoff, n, which 
determines the level of noise to be absorbed, and the parameter Err for esti- 
mating the precision. 

1. Let us take a trial set of amplitudes B n (n = . . . N — 1), that is the real 
and the imaginary parts Re B n , Im B n and the absolute values \B n \. 

2. Consider the evolution of B n described by equation ([h]) for the initial 
values given by the set B n , that is B n (t = 0) = B n (n = 0. . . N — 1) for 
the time interval At, and take the set of amplitudes B n = B n (t = At). 

3. Consider the set of absolute values \B n \, and find the index M for the 
maximum value, Bm, of \B n \ (n — . . . N — 1). 

4. Find the indices L and R such that \Bl\ and \Br\ are both local minima 
of the set \B n \, and the constraints 

(i) \B L \ <IIand \B R \ <n, 

(ii) L < M < R, 

(iii) L, R are closest to the index M, 
are verified. 

5. Consider the new set of amplitudes B n (t) 




otherwise 



L < n < R 



where 



2 



A= ^2\B n (At) 



_n—L 
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6. Consider B n {At + T) which are the solution to equation (|T(]) for the initial 
values given by B n (At) and the time interval of integration equal to T. 
Verify that the mean square deviation 

^2[\B n (At + T) -B n {At)} 2 , At<T 
n 

is less than the accuracy level Err accepted. In case it is not met, set 

B n = B n (t = 0) 

else exit. 

7. Return to step 2. 



In case of solitons with large support, or long wavelengths, much larger than 
the spacing a, we can use the usual approach and look for solitons of the form 
@, that is 

Y(x, t) = e l ^ na ~^ ^(na - vt) 

On substituting B n (t) given above in equation (|l0|), we obtain the two standard 
equations for the imaginary and the real parts of B n (t) . In the long wavelengths 
limit, the imaginary one results in the equation for the velocity of soliton, which 
in our notation reads: 

v = 2a ■ up ■ sin(a • q). 

The real part gives a nonlinear functional equation for the amplitude tp, which 
can be solved by the familiar Newton method. 

It is important that, for long wavelengths, the RF-algorithm brings about 
the same results as the Newton method indicated above. 
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Figure Captions 



Fig.l 

Schematic structure of the double helix of DNA. 
Fig.2 

Amino/imino and keto/enol forms of purine and pyrimidine. 
Fig.3 

Thymine- Adenine and Cytosine- Guanine pairing. 
Fig.4 

Pairing of cytosine, normal form, and adenine, imino form. 
Fig.5 

Pairing of thymine, enol form, and adenine, imino form. 
Fig.6 

Uracil-uracil pairing. 
Fig.7 

Pairing of adenine- (brom-uracil) and guanine- (brom-uracil). 
Fig.8 

Typical moving soliton for W = 0.75, A = 0.001, velocity 1325.28 base pairs 
per msec, the period of time spent 0.533 msec. The solid line shows \A n \ 2 , and 
the thin lines indicate the real and the imaginary part of the amplitude A n . 

Fig.9 

Breather, or still soliton. Typical moving soliton for W = 10, A = 0.5, velocity 
base pairs per msec, the period of time spent 0.501 msec. The solid line 
shows \A n \ 2 , and the thin lines indicate the real and the imaginary part of the 
amplitude A n . 

Fig.10 

Radiation emitted from the motionless central peak during the period of time 
0.018, for W — 5 and A = 0.5; the velocity of side waves 1833.3 base pairs/msec, 
distance traveled 33 base pairs. The solid line shows \A n \ 2 , and the thin lines 
indicate the real and the imaginary part of the amplitude A n . 

Fig.ll 

Moving soliton for W — 0.75, A = 0.075, velocity 1335.46 base pairs per msec, 
the period of time spent 0.530 msec, distance travelled 707 base pairs. The solid 
line shows |A„| 2 , and the thin lines indicate the real and the imaginary part of 
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the amplitude A n . 
Fig.12 

Moving soliton for W = 2, A = 0.1, velocity 851.2 base pairs per msec, the 
period of time spent 0.577 msec, distance travelled 491 base pairs. The solid 
line shows \A n \ 2 , and the thin lines indicate the real and the imaginary part of 
the amplitude A n . The values of W, A are close to the borderline (see Fig. 17), 
dividing the region of stable solitons from the unstable ones. 

Fig.13 

Partial self-focusing of an initial low amplitude distribution on a peak for a 
period of time 0.002 msec, for W = 1, A = 0.5. 

Fig.14 

Pairing of thymine and 2'-Deoxyisoguanosine. 
Fig.15 

Filter algorithm for finding solitons. The solid line indicates the part of the 
excitation to be preserved, and the thin line the does the part to be cut off, the 
dashed line means the precision. 

Fig.16 (a)-(d) 

Sets of W and v, which allow for soliton solutions, at fixed A: 

(a) for A = 

(b) for A = 0.1 

(c) for A = 0.15 

(d) for A = 0.2 

The solid line shows still solitons, or breathers, and the shaded area indicates 
moving solitons. We take into account only solitons of size less than 100 base 
pairs. 
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